function [H,row_idx] = smooth_C1(dof_map,V,T,TE,ET,d)
% function [H,row_idx] = smooth_C1(dof_map,V,T,TE,ET,d)
% get the global smoothness conditions 
idx_mat = zeros(3,d,3); % the column idx 1 2 3 -1 -2 -3
idx_c0 = zeros(d+1,3);
cr = cr_pattern(d);
for i = 1:3
    idx_c0(:,i) = cr_indices(0,d,i,cr);
    idx_mat(i,:,1) = cr_indices(1,d,i,cr);
    line = cr_indices(0,d,i,cr);
    idx_mat(i,:,2) = line(1:d);
    idx_mat(i,:,3) = line(2:d+1);
end

E_in = find(ET(:,2)~=0); nE_in = length(E_in);
%row_idx = zeros(nE_in*(2*d+1),1);
row_idx = zeros(nE_in*d,1);
%i = zeros(nE_in*(4*d+1),1); j = zeros(nE_in*(4*d+1),1); s = zeros(nE_in*(4*d+1),1);
i = zeros(nE_in*3*d,1); j = zeros(nE_in*3*d,1); s = zeros(nE_in*3*d,1);
pos = 1; row = 1;
v_idx = [1 2 3;2 1 3;3 1 2];
for k = 1:nE_in
    eg = E_in(k);
    tri = ET(eg,1); nei = ET(eg,2); % find two triangles who have eg.
    eg2tri = find(TE(tri,:)==eg);  % local index of eg in tri
    eg2nei = find(TE(nei,:)==eg); % local index of nei in tri

    % first set c0 condition:
%    i(pos:pos+d) = row:row+d;
%    j(pos:pos+d) = dof_map(idx_c0(:,eg2tri),tri);
%    s(pos:pos+d) = 1; pos = pos + d + 1;

    % second set c1 condition:comput the barycentric of v4(b1,b2,b3) relative to tri
    v1 = v_idx(eg2tri,1);    v2 = v_idx(eg2tri,2);    v3 = v_idx(eg2tri,3);
    A = [1 1 1;V(T(tri,v1),1) V(T(tri,v2),1) V(T(tri,v3),1);V(T(tri,v1),2) V(T(tri,v2),2) V(T(tri,v3),2)];
    b = A\[1;V(T(nei,eg2nei),1);V(T(nei,eg2nei),2)];
    % for each coordinates set c1 condition
    for l = 1:3 % the idx is always from v2->v3
%        i(pos:pos+d-1) = row+d+1:row+2*d;
        i(pos:pos+d-1) = row:row+d-1;
        j(pos:pos+d-1) = dof_map(idx_mat(eg2tri,:,l),tri);
        s(pos:pos+d-1) = b(l);
        pos = pos + d;
    end
    
    % both c^0 and c^1
%    if mod(eg2tri+eg2nei,2)==1 % the same direction with v2->v3
%        row_idx(row:row+d) = dof_map(idx_c0(:,eg2nei),nei);
%        row_idx(row+d+1:row+2*d) = dof_map(idx_mat(eg2nei,:,1),nei);
%    else  % the different direction in neighbour
%        row_idx(row+d:-1:row) = dof_map(idx_c0(:,eg2nei),nei);
%        row_idx(row+2*d:-1:row+d+1) = dof_map(idx_mat(eg2nei,:,1),nei);
%    end
%    row = row + 2*d + 1;

    % only c^1
    if mod(eg2tri+eg2nei,2)==1 % the same direction with v2->v3
        row_idx(row:row+d-1) = dof_map(idx_mat(eg2nei,:,1),nei);
    else  % the different direction 
        row_idx(row+d-1:-1:row) = dof_map(idx_mat(eg2nei,:,1),nei);
    end
    row = row + d;   
end
dim = max(max(dof_map))-min(min(dof_map))+1; pos = pos-1; row = row-1;
H = sparse(i(1:pos),j(1:pos),s(1:pos),row,dim) - sparse((1:row)',row_idx,ones(row,1),row,dim); 

function cr = cr_pattern(d)
% the matrix store the indices of degree d
% much alike the function asce_pattern and desc_pattern
% the local dof index is stored in the up-right triangle.
rows = d+1;
cols = rows;
cr = zeros(rows,cols);
count = 1;
for line = 1:rows
    for row = line:-1:1
        cr(row,line-row+1) = count;
        count = count + 1;
    end
end

function line_dof = cr_indices(r,d,idx,pattern)
% to get the r's line of dof_map of degree d parallel to edge i
% it is specially useful when treating smoothness condition and boundary
% condition! 
% the output parameter mat is the matrix P in my notes.
switch idx
    case 1
        line_dof = (((d+1-r)*(d-r)/2+1):((d+2-r)*(d+1-r)/2))';
    case -1
        line_dof = (((d+2-r)*(d+1-r)/2):-1:((d+1-r)*(d-r)/2+1))';
    case 2
        line_dof = pattern(r+1,1:d-r+1)';
    case -2
        line_dof = pattern(r+1,d-r+1:-1:1)';
    case 3
        line_dof = pattern(1:d-r+1,r+1);
    case -3
        line_dof = pattern(d-r+1:-1:1,r+1);
    otherwise
        line_dof = [];
end
%then put the corresponding value to a matrix;
% i = (1:d-r+1)';j = line_dof; s = ones(d-r+1,1);
% mat = sparse(i,j,s,d-r+1,(d+1)*(d+2)/2);